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I. INTRODUCTION 

Reacting complex flows are ubiquitous in the atmosphere, the oceans, and the interior 
of our planet: Stratospheric ozone chemistry is perhaps the most well known example of 
chemical dynamics occurring in a stirred geophysical fluid. But the fate of many other 
reacting tracers is also influenced by mesoscale eddies in the ocean, or by convection in the 
Earth mantle, to name just a few cases. Biological population dynamics is also a kind of 
"chemical reaction". The interactions of nutrients and plankton species in the sea are the 
first steps in the ocean food chain, and they are among the most important ingredients for 
the understanding the CO2 interchange between atmosphere and ocean. 

Tracers stirred by fluid motion are known to develop strong inhomogeneities, usually 
in the form of filamental features, arising from a kind of variance cascade from the forcing 
scale towards smaller scales. These structures are poorly resolved in global atmospheric and, 
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especially, in oceanic models. They provide however sensitive mixing mechanisms and they 
are in some sense "catalists" |I[ that enhance the chemical or biological activity occurring 
in geophysical flows. We mention as examples the increase of ozone depletion produced by 
the filamental structure of chlorine filaments [Q], or the increase in biological production in 
small-scale ocean structures [[J. 

In this set of Lectures, we will present some theoretical ideas for the study of structures 
arising in advected reacting tracers. The emphasis will be in the modifications that chemical 
or biological activity introduces in reactive patterns with respect to the ones appearing in 
passively advected substances. 

Examples from plankton populations dynamics will be used through these Lectures. It 
should be mentioned, however, that if one looks into the literature, the diversity of observa- 
tions related to plankton distributions is so great that some authors, rather pessimistically, 
argue that "... One should expect no general results that describe in all (or most) cases how 
'the biology' modifies a spatial pattern that has arisen solely from the dispersal of species 
..." ||. We use here the biological experimental observations more as a motivation and 
illustration, than as a detailed check of a particular theory. Our aim is to point out simple 
and robust mechanisms that can be useful in classifying the large number of observations 
into a more reduced set of different scenarios. Later quantitative test would need synoptic 
measurements of the biological variables and of the hydrodynamic flows. This would be 
certainly more accurately done in laboratory chemical experiments than in open sea biol- 
ogy. Most of the results presented here are of direct application to chemical reactions of the 
decaying or of the excitable type. Nevertheless, we mention that comparison of some of our 
theoretical ideas with simultaneous satellite data of both a reactive tracer (sea temperature 
in interaction with air temperature) and hydrodynamic flow has been very recently initiated 

In Section [TI] we will briefly review some aspects of plankton dynamics that will be useful 
in the following. The paradigm of chaotic advection, i.e. the chaotic transport of fluid 
parcels by smooth large-scale flows, will be the one more frequently used here to discuss 
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flow effects, and it is described in Section |T|. Sections [TV] and [V] describe results obtained 
for two rather general classes of chemical or biological reactions in chaotic flows: the first- 
order reaction, or linearly decaying tracer, and reactions of the excitable type. Finally (Sect. 
|VID , some warnings will be given about the difficulties in modeling discrete individuals (such 
as planktonic organisms) in terms of continuous concentration fields. 



II. PLANKTON, PATCHES, AND BLOOMS 
A. Background 

Plankton is the generic name given to a huge variety of aquatic organisms, both from 
marine and from freshwater environments, comprising from viruses with less than 200 nm 
in size, to crustaceans or even fish larvae, with sizes above the centimeter. What they have 
in common, and distinguishes them from other major group called nekton, is their inability 
to overcome the major ambient hydrodynamic currents. In this sense, they are strongly 
influenced by the existent hydrodynamic weather ||. It is often assumed that they are 
passively transported by the marine currents. This is probably true for the smallest species, 
but it should be said that most of these organisms present flagellae and other natatory 
organs that give to them some mobility, especially in the vertical direction. In addition to 
interactions with the hydrodynamic environment different plankton species interact among 
them in a variety of ways, such as predator-prey interaction, parasite-host relationship, 
competition for resources, etc. The biological interaction between species and with the 
nutrient substances dissolved in water is formally equivalent to a chemical reaction dynamics. 

Any attempt to overview here some major theme in plankton biology and its interaction 
with hydrodynamics would result very partial and incomplete, and we refer the interested 
reader to some monographs that turn out to be very readable. We just mention that a 
major division between plankton species arises between the ones that are able to assimilate 
inorganic material by photosynthesis, the phytoplankton, and the species that necessarily 
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graze on the above, the zooplankton. In the following we will briefly describe just two 
phenomena, plankton patchiness and plankton blooms, extremely stimulating from the per- 
spective of cross-disciplinary research, and for which some ideas from Nonlinear Science may 
be of relevance. 

Plankton patchiness |J is the term used to label the large inhomogeneity observed in 
plankton distributions. Earlier observations from land and from ships reported water masses, 
of sizes around say 10 km, distinctly colored by the presence of different species of planktonic 
organisms. The search for biological or physical mechanisms singling out a characteristic 
patch size, despite mixing with the surrounding water stimulated the first theoretical ap- 
proaches (Ty|ll| • Modern satellite observations (T^j reveal that the range of sizes of plankton 



structures is much broader. The largest features, thousand of kilometers wide, are associated 
with major ocean gyres and currents. Upwellings of nutrient-rich water from deeper ocean 
layers, and river discharges, fertilize the upper layer of the ocean and induce associated 
plankton structures. In the mesoscale range, fronts and eddies are seen as major players, 
both in bringing to the surface nutrients from lower layers, and in providing horizontal 
transport and stirring processes. 

A closer look at the inhomogeneities in plankton distributions can be performed by ana- 
lyzing water samples along a ship trajectory (or transect). Higher spatial resolution can then 
be attained and the result is an extremely intermittent distribution (Fig. [I], lower panel). 
At this point, it is pertinent to ask about the relevance of the fact that plankton is a liv- 
ing substance, or if perhaps the same kind of distributions can be found for non-biological 
tracers. A partial answer to this question is given by Fig. [I], in which distributions of phys- 
ically controlled tracers (temperature and salinity) are compared with the phytoplankton 
concentration, indicated by water induced fluorescence, at a resolution of about half a kilo- 
meter. Although correlations can be appreciated, evident differences are seen between the 
distributions of both kinds of tracers. A brief walk on the literature shows that the distribu- 
tions may look very different form place to place in the ocean, and that their characteristics 



may also change in time In general one can say that they will be affected by both 



the characteristics of the fluid flow that transports them, and by the chemical or biological 
interactions in which they are involved. 



Latitude II 8 Nov 1995 



U3 
CO 

'2 



25.30 
25.20 
25.10 
25.00 

37.20 
37.15 
37.10 
37.05 
37.00 

0.014 
0.012 



o 0.010 

w 
a, 
s-. 
o 

p I 







-J 






— 1 1 1 1 1 1 1 








1 





100 



200 

Distance (km) 



300 



FIG. 1. Temperature (in degrees C), salinity (in Practical Salinity Units) and induced fluorescence (in 
arbitrary units) measured at 5 m depth along the path of the research vessel BIO Hesperides, across the 
southern Atlantic tropical gyre during the Latitude II campaign |l5| . Fluorescence identifies chlorophyll and 
is thus a proxy for phytoplankton concentration. The starting point of the transect shown was at 31.65 W, 
16.48 S, and the final point at 32.83 W, 19.43 S. 

The irregular nature of the observed distributions immediately suggest characterization 
in terms of scaling exponents. This has traditionally been done in terms of the power spec- 
trum of concentration fluctuations ||14|| , since a comparison with Kolmogorov-like turbulence 
theories was then direct. In situations in which the phytoplankton concentration spectrum 
displayed a large and well defined power-law regime (which is not the case of the data in 
Fig. |l]) it decayed as k~@, with k the wavenumber. The scaling exponent (3 results to be 
larger than one. We prefer to discuss the scaling properties of concentration fluctuations 
in terms of structure functions, since a more complete characterization of intermittency 



properties can be done. These quantities, and their relation with the power spectrum, are 
introduced in Sect. [TV]. 

A second kind of inhomogeneity found in plankton distributions is temporal inhomo- 
geneity associated to plankton blooms. Seasonal blooms are episodes of rapid phytoplank- 
ton growth that occur in mid-latitude seas at the beginning of spring. This is followed by 
zooplankton (and also bacterial) growth that consumes the available phytoplankton until 
concentrations return back to the initial low level. The whole cycle finishes with the summer. 
Sometimes a smaller version of the bloom repeats in autumn. The mechanism for seasonal 
blooms seems to be understood, at least qualitatively: the beginning of ocean stratifica- 
tion in spring reduces the depth of the layer in which phytoplankton is mixed by turbulent 
motion, the mixed layer, so that these photosynthesis-capable organisms spend more time 
exposed to solar light. Fast growth then occurs. The subsequent growth of predators, and 
the consumption of nutrients in the upper layer, ends the bloom. Destratification, with 
the destruction of the seasonal thermocline in autumn, brings again to the surface some 
nutrients from the deeper layers, so that phytoplankton activity may increase again. 

There is another kind of blooms that occur more localized in space, and not with a 
seasonal periodicity, but occasionally and associated to local warming or local increase of 
nutrients in the water. They are called sometimes red tides since the large amount of 
organisms in water may give to it some distinct coloration. Time scales for the initial 
development of the bloom are in the range of days, and the end usually arrives in a time of 
the order of weeks or month. 

B. Models 

Modeling the above spatial and temporal evolutions of planktonic populations is a chal- 
lenge since many physical and biological processes are involved, some of them poorly un- 
derstood. A variety of levels of description are available, from detailed models based on 
individual organism behavior, to mean field approaches that consider the whole sea as an 
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homogeneous soup without spatial features. Another aspect of modeling that should be 
addressed is the degree of aggregation at which groups of species are represented by the 
same collective variables. It should be said that these starting points, the election of the 
adequate levels of description, and the relationship between different degrees of detail, is an 
important subject by itself, and its complexity is just beginning to be grasped [|T^,|T^]. More 
on this in Sect. IVl 
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FIG. 2. The processes in the NPZ models. The letter in parenthesis identifies the process with the term 
representing it in Eqs. (|]). 



A class of models simple enough to allow some explicit analysis, but still retaining some of 
the complexities of real-life trophic interactions is the one known as the NPZ class of models 
(from Nutrients, Phytoplankton and Zooplankton). We describe it here for illustration, and 
because it will be of use in some of the results presented later. The processes included 
within this framework are displayed in Fig. |2], and its mathematical modeling is as follows 
[|T8| , [T9[1 : in each portion of well mixed fluid, the time evolution of the amount of nutrients 
(N), phytoplankton (P), and zooplankton (Z) is ruled by 

^ = F N (N, P, Z) = s(N — N) — f(N)P + cP + ag{P)Z + eZ 
at 



F P (N, P, Z) = f(N)P -qP- g(P)Z 



dP 
~dt 

— = F Z (N, P, Z) = kg(P)Z - mZ* 
at 



(1) 



The nutrient input Nq affects the nutrient concentration in the mixed region, and this 
influence is transferred to phytoplankton by the second equation, and to zooplankton by the 



third. Finally, zooplankton either dies or is consumed by higher organisms, which leads to 
the last term in the third equation. The exponent r in this last term is usually taken to be 
either r = 1 or 2, depending on the dominant kind of predation on zooplankton one wants 
to model. 

For the description of the spatial dependence, the framework of Advection- Reaction- 
Diffusion equations is frequently used. They are partial differential equations in which 
the organisms and their nutrients are described by continuous distributions subjected to 
advection by a velocity field v(x, t), and to diffusion, in addition to the interactions ([[]). 
The general form of the equations is 

^ + V-(vQ) = F 4 (iV,P,Z) + AV 2 Q . (2) 

The index i takes three values, so that Fi is either Fn, Fp, or Fz, and the variables Cj are 
either N, P, or Z (and now these symbols denote concentrations per unit volume). If the 
flow is incompressible (V • v = 0) then the second term in the left- hand- side can be written 
as v • VCj. At scales large enough, so that the effects of diffusion can be neglected, Eq. (0) 
admits a Lagrangian representation that in the incompressible-flow case coincides with the 
system in Eq. ([TJ), where now N, P, and Z are either the concentration or the amount of 
substance contained in a fluid element that moves with the flow velocity. The trajectory 
r(t) of this fluid parcel thus satisfies 

-=v(r, i ) . (3) 

It may seem that Eqs. (P and (|3|) represent uncoupled dynamical systems. This would 
be the case if none of the parameters in ([]]) is spatially dependent. But if one of them 
varies from point to point, it should be evaluated at the position of the fluid element at 
time t, that is, at r(t). This couples the chemical and the transport dynamical systems. For 
example, the parameter Nq may represent an inhomogeneous nutrient input iVo(x) into the 
upper ocean mixed layer by a localized upwelling, and thus the fluid element following the 
trajectory r(t) would experience a time-dependent input AT (x = r(t)). 
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The model is still largely undefined until one specifies the parameters and, particularly, 
the response functions / and g, which really contain information about the species interac- 
tion. A rather common form for the function /, experimentally established by Monod, is 
the so-called Michaelis-Menten form 

which may be also used for the grazing function {g{P) = aP/(P e + P)). This form may be 
justified by arguments taken from enzyme kinetics ||20|| . Any function that, like (£|), presents 
a linear increase close to the origin followed by saturation at large N receives the name 
of Hollings type-II response function. Another commonly used Hollings type-II response 
function is the Ivlev one: g(P) = a(l — e xp ). Hollings type-Ill response functions differ 
from the above by its behavior at small argument: the linear increase is replaced by a slower 
one, indicating that predators will not be interested in a too small prey concentration. An 
example of this type is 

Large uncertainties exist in the form of the response functions and in the numerical 
values of the parameters involved in ([!]) for particular ecosystems. We mention however 



that quantitative determination can be achieved under laboratory conditions [gj] . 

Many variations and simplifications of model (jl[) can be found in the literature. It 
is rather common to reduce the modeling to the predator-prey competition of zooplankton 
and phytoplankton. In this case the limitation imposed on the phytoplankton growth by the 
nutrient concentration is included by replacing the second equation in (U) by one containing 
a logistic term: 

F P (P,Z)=aP(l-^)+... (6) 



C 

The parameter C is called the carrying capacity, and represents the maximum amount of 
phytoplankton the fluid parcel can support. 
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C. Some proposed mechanisms for patchiness 



Probably, among the first theoretical approaches used to address the formation of plank- 



ton patches is the one contained in papers by Skellam [10], and by Kier stead and Slobodkin 
1 1| (the so-called KISS theory). They considered the free growth (at rate fi) of phyto- 
plankton (with concentration P(x, £)), in a finite region of water suitable for growth, with 
plankton escaping the suitable zone by diffusion: 

BP 

^- = ^P + DW 2 P . (7) 

This is a very particular case of (§). It is assumed that P = outside the region. If the 
patch is too small, the growth can not overcome the diffusive escape flux, and the patch dies. 
This can be seen, for example, via the exact one-dimensional solution of fl7|), with boundary 
conditions such that phytoplankton disappears outside the region [0, L] of suitable conditions 
for growth (i.e. P(x = 0,t) = P(x = L,t) = 0): 

P(x,t) = ±A n sin(^)e^~ D ^ 2 ) . (8) 



71=1 



The constants A n are fixed by the initial conditions. It is clear that if L < L c = iryD/fM no 
growth occurs and the patch dies, so that L c can be identified with the minimum patch-size 
able to support growth. 

There is a number of criticisms than can be posed to this approach. The first one is 
that the model is incomplete, since once growth begins it continues without limit. Nonlin- 
ear saturation and interactions with predators would be needed to stop this. The diffusion 
coefficient D is certainly not the originated from the Brownian motion of the organic parti- 
cles, since this would be irrelevant to processes above, say, the millimeter scale. It is rather 
a turbulent eddy-diffusion coefficient aimed to represent in an averaged way the effect of 
dispersion by the turbulent flow, with values that would depend on the observation scale as 
found for example in |2!| . One should ask, thus, how can the plankton disperse in such way, 
whereas the nutrients remain static within the fixed parcel, despite they live in the same 
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flow. In any case, if one introduces a growth rate for phytoplankton of the order of one day 
(/i « 10 -5 s -1 ) and an eddy diffusivity appropriate for the 10-100 km range (D « 100 m 2 /s) 
P2|| , one finds L c « 10 km, which seemed to compare well with the earliest observations of 
plankton patches. 

An interesting development that has some common ideas with the KISS approach, but 



introduces a new and important element, is described in p3| . The idea is that what will con- 
trol dispersion of a patch is not only eddy diffusivity, but also the geometric characteristics 
of the mean flow. It turns out that if an incompressible fluid flow induces dispersion in one 
direction it necessarily produces concentration in another, to conserve the fluid volume. The 
situation is especially clear in a purely two-dimensional flow, where the rate of expansion 
of a small fluid blob in one direction should be exactly the same as the rate of contraction 
in another. Thus, an initially circular patch will take the form of a filament locally aligned 
with particular directions in the flow, the unstable manifolds of the flow hyperbolic points. 
These manifolds have been nicely observed in laboratory experiments |24j] , and phytoplank- 



ton filaments on the sea surface have been observed from space Reference |33| proposes 
the following model to describe the transverse profile of a phytoplankton filament: 

_-A,_ =„(*),> + ,>_ . (9) 

There are two differences with the KISS model. One is the possible time-dependence of 
the growth rate n(t), that is a simplified linear way of modeling nonlinear interactions and 
predation on phytoplankton after the first growth stages. The other is the advective term 
—Xxd x P, associated with a velocity field v x (x) = —Ax that models a local strain that, in the 
direction perpendicular to the filament, compresses it. A solution to Eq. (|9]) that is stable 
and attractive can be found analytically ||23|| : a Gaussian with time-dependent height but 



fixed width w = yD/X. Thus, the lateral scale of the filament is not controlled by some 
externally imposed size of the suitable water, but by the competition between diffusion 
and advection. It is somehow surprising that the biological growth rate does not affect the 
filament width in this model. As we will see in Section [V] this feature is not shared by other 
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models. 

A different mechanism that has been invoked as possible source of plankton patchiness 
|26| is the Turing mechanism. It appears generically when there is competition between 
a self-reproducing prey and a predator that diffuses faster than the prey, and leads to 
patterns with a characteristic periodicity J20J. Although diffusion coefficients of the organic 
particles have roughly the same very small value for both predator (zooplankton) and prey 
(phytoplankton), the former has usually a larger vertical mobility. If there is some vertical 
shear in the water column, the organisms with larger vertical mobility will experience a 
more variable velocity field, and thus they will be subjected to larger dispersion. If one 
models such enhanced dispersion by an effective diffusion coefficient, it would be larger for 
zooplankton than for phytoplankton, and then the Turing mechanism can be at work. The 
analysis in |27]] implies that for models of the type (||), the Turing mechanism does not occur 
in the relevant range of parameters. But this may not be the case for modifications of this 
class |[28|| . In addition, related mechanisms, that would lead also to plankton distributions 
spatially periodic, may appear in the presence of differential motion between predator and 
prey [f£|. 

Horizontal stirring is another mechanism that produces strong inhomogeneity PUI . This 
is the one that will be addressed in the following. According to , the attractors of models 
of the form (]l|) are, for most of the relevant parameter range, fixed points representing stable 
species coexistence. In most cases, at least if type III Hollings functions are used, excitable 



behavior f20[ occurs close to these equilibrium points. This motivates us to consider the 
effect of stirring on population models displaying either simple relaxation to fixed points 
(Section [IVD or excitable behavior (Section 0). Population oscillations also occur at realistic 
parameter values ||18| , |2~Ij| , but we will not consider this situation here. Before going to the 
results for the chemical/biological dynamics, we summarize in the next section some results 
about chaotic advection, which is the kind of fluid flow considered in the developments 
presented here. 
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III. THE PARADIGM OF CHAOTIC ADVECTION 



When dealing with transport processes in complex fluid flows, the concept of turbulence, 
the unsteady and irregular kind of flow in which motion occurs in a large range of scales, 
is the paradigm that most often comes to mind. Fluid elements, and transported particles, 
follow intricate trajectories in such velocity fields. It was recognized some time ago pTT 



that a turbulent velocity field is not necessary to produce chaotic fluid-element trajectories. 
In three spatial dimensions the dynamical system (^) has chaotic trajectories P2| , |3"3"|| even 
for very simple steady and smooth velocity fields v(x). In two-dimensional flows, a simple 
periodic time dependence in v(x, t) is enough to induce generically chaotic trajectories in 
This situation of chaotic trajectories produced by a simple and regular velocity field 
is called chaotic advection or Lagrangian turbulence. The essential characteristic of chaos 
is sensibility to initial conditions. This means that fluid elements initially very close would 
follow diverging trajectories. The exponential rate of separation is given by the flow Lya- 
punov exponent A. In an incompressible fluid, contraction must occur in another direction 
to conserve volume. The clearest situation is in hyperbolic two-dimensional flows, where one 
can define at each point an expanding direction and a contracting direction, with both the 
contraction and expansion rate occurring at the same exponential rate A. As a consequence, 
portions of fluid initially compact become stretched in long and thin filaments, and are also 
repeatedly folded to remain inside the system. As a result, chaotic advection produces a 
"cascade" of inhomogeneities from large scales to smaller and smaller scales, in the same 
qualitative way as a fully turbulent velocity field will do. Quantitative details are, however, 
different ||34|| . It turns out that the patterns produced by advection in smooth velocity fields 
are more singular than the ones produced by singular velocity fields such as Kolmogorov 
turbulence. 

Experimental realizations of (nonturbulent) chaotic advection can be performed in small 
containers with high viscosity fluids, so that Reynolds number remains small and the flow 
laminar. It is now recognized that the regime of chaotic advection also applies to turbulent 
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flow at scales smaller than the Kolmogorov scale ( 1% = (^/e) 1 / 4 , where v is the kinematic 



viscosity and e the energy dissipation rate). This is the so-called Batchelor regime [35 



The Lyapunov exponent in that range is of the order of the mean strain A ~ \j e l v - F° r 
water [y « 10 -6 m 2 s _1 ) at conditions typical in the upper mixed layers of oceans and lakes 
(e ~ 10~ 7 W kg -1 , [[3(J), Ik is between the millimeter and the centimeter. Thus the Batchelor 
regime is certainly relevant for plankton processes such as feeding, aggregation, encounter, 
etc. It is however a range much smaller than the structures shown in Fig. [L] or the ones 
visible from satellites. A third situation where chaotic advection has been used with success 



is in the modeling of the largest scales of geophysical flows |37] . Smaller motions would then 
have to be represented by some eddy diffusion. 

In the following, we focus on the situation of chaotic advection. This would allow us to 
progress in our objective of identifying mechanisms that can contribute to classify obser- 
vations in robust scenarios. Chaotic advection is a framework simple enough to obtain a 
number of explicit results that turn out to be rather insensible to model details. Of course, 
explanation of particular geophysical observations would require additional discussion on 
whether or not the paradigm of chaotic advection applies to the scales considered. 



IV. THE DECAYING TRANSPORTED SUBSTANCE: A MODEL FOR STABLE 

DYNAMICS 

Most of the simpler models of chemical dynamics, and models such as ([!]) in a large 
range of parameters, evolve in such a way that, at long times, the different species reach an 
equilibrium at which they coexist at some fixed-point concentration values. When some of 
the coefficients in Eq. ([!]) are functions of space they act (when coupled with the transport 
equation (|3])) as time-dependent source or sink forcing terms that disturb such equilibrium. 
Even in this forced case it may happen that the concentration values at each fluid element 
tend to relax to a unique time-history, determined by the fluid-element trajectory. When 
this happens we say that we have a stable chemical dynamics. In a sense, the concentra- 
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tions try to relax to the local-equilibrium values that would correspond to the values of 
the parameters found along the fluid-element trajectory. The mathematical way to check 
this stability is by initializing a particular fluid element with two slightly different sets of 
concentration values, and monitor how these values evolve under the coupled chemical and 
transport dynamics (|I|) and (|3|). The exponential rate of divergence between these two dif- 
ferent initial concentrations, under the same fluid trajectory, defines the so-called chemical 
Lyapunov exponent, Ac- If this quantity is negative, it identifies convergence to a common 
evolution, and this is the case of stable chemical dynamics. The usual ergodic properties 
would guarantee that the same value of the chemical Lyapunov exponent will be obtained 
for almost all fluid trajectories and initial concentrations. 

The simplest case of stable chemical dynamics is the linearly decaying passive scalar, or 
first order reaction. It consists in the evolution of a single concentration C(x, t) under a 
linear decay at rate b: 

BC 

— + v(x,t) ■ VC = 5(x) -bC + D\7 2 C, (10) 

A space-dependent forcing consisting in a source S(x) of the substance is included, to main- 
tain a non-trivial concentration field at long times. 

It is easy to check (by using the Lagrangian representation valid for D — > 0) that 
Ac = —b. The simple dynamics (0) can be considered either as an approximation to more 
complex chemical or biological evolutions with stable chemical dynamics, or as a description 
of simple specific processes such as spontaneous decomposition of unstable radicals, evolu- 
tion of a radioactively or photochemically decaying substance, or relaxation of sea-surface 
temperature towards atmospheric values ||. It also describes the concentration of the C 
reactive in the binary reaction B + C — > D when the B concentration is maintained constant. 

The first systematic study of model (|T0"| ) seems to be in ]58 |. The author analyzed, by 



a spectral element method, the power spectrum in a number of threedimensional turbulent 
regimes. In the Batchelor regime a power-law is found with an exponent depending on 
the quotient between the decay rate b and the mean strain rate, which is equivalent to 
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a flow Lyapunov exponent A. Here we will obtain an essentially equivalent result within 



the framework of chaotic advection p9| , ^0| . In addition, intermittency corrections and the 
differences between open and closed flows will be discussed. Our results are first obtained 
in terms of the behavior of the concentration differences between spatially close points 
(5C = C(x + <5x, t) — C(x, t)) and later in terms of the structure functions, that are the 
moments of SC averaged in space. 



Exact results in |41j for the decaying scalar in a closed random Kraichnan flow in the 
Batchelor regime confirm the generality of the mechanism found: chaotic advection and 
decaying chemistry produce singular patterns that can be characterized by scaling exponents 
depending on the chemical decay rate and on the flow Lyapunov exponent (more precisely, 
on the finite-time Lyapunov exponent distribution). We will present here also results for 
nonlinear plankton models, showing that the results of the linearly decaying chemistry can 
be translated to the case of nonlinear stable chemistry at the smaller scales simply by 
replacing the decay rate b by the absolute value of the chemical Lyapunov exponent |Ac|. 
Detailed justification of this can be found in H2j for multiple species in the framework of 
chaotic advection, or in [I3| for a single concentration in the framework of the random 
Kraichnan flow. It is important to mention that in the present framework the relevant time 
scale affecting the distributions structure is identified to be Ac, and not other time scales 
frequently used in the literature such as the phytoplankton linear growth rate. 



A. The smooth-filamental transition for the decaying passive scalar in a closed flow 

We consider (|T0| ) with a two-dimensional, incompressible, smooth, and non-turbulent 
velocity field. Chaotic advection is obtained generically if a simple time-dependence is in- 
cluded in v(x, t). We assume that diffusion is weak and transport is dominated by advection. 
Thus one expects that the distribution on scales larger than a certain diffusive scale is not 
affected by diffusion, and set D = 0. In this limit the above problem can be described in a 
Lagrangian picture by the pair of dynamical systems ([|) and 
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dC 

— = S[x = r(t)]-bC, (11) 

where the solution of (§) gives the trajectory of a fluid parcel, r(t), while (|TT|) describes the 
Lagrangian chemical dynamics in this fluid element: C(t) = C(x = r(t),t). 

To obtain the value of the chemical field at a selected point x at time i one needs to 
know the previous history of this fluid element, that is the trajectory r(t) (0 < t < t) with 
the final condition r(t) = x. This can be obtained by the integration of (|3]) backwards in 
time. Once r(t) has been obtained, the solution of ( |TTD is 

Che, i) = C[r(0), 0}e- M + f S[v(t)\e~ h ^ dt. (12) 

Jo 

Then one can obtain the difference at time i of the values of the chemical field at two 
different points x and x + 5x separated by a small distance <5x in terms of the difference 
8C[r(t),t\ 8r(t)\ = C[r(t) + 8r(t),t] - C[r(t),t] for < t < i, namely: 

8C(r,t;8x) = 6CM0), 0; <5r(0)] e - fe * + /* 8S[r(t); 8r(t)]e- b ^dt (13) 

Jo 

where 8r(t) (0 < t < t) is the time-dependent distance between the two trajectories that 
end at x and x + <5x at time i, and 8S is the difference of the source term at points r(t) and 
r(i) + 5r(t). In the following we omit the first term in the right-hand-side of (0) since it 
always disappears at long times. 

If the trajectory is chaotic, we have in the backwards dynamics, for t < and large, 
|<5r(t)| ~ |5r(0)|e _A *, where A is the positive value of the Lyapunov exponent along the 
trajectory. In flT3|), 8r is obtained backwards starting from 5x at t = i. In this case: 

8r(t) w 5x e A(f -* ) , t < i . (14) 

We note that there is a particular direction for the orientation of the initial displacement 5x 
along which the trajectories converge instead of diverging (this is the contracting direction 
in the backwards flow which is the expanding one in the forward dynamics) . For <5x oriented 



along it, —A in ([14]) should be replaced by A. For hyperbolic dynamical systems, the value 

17 



n ■ VCYr, t) » n ■ / V5[r(t)]e (A ~ 6)(f -* ) cit . (16) 

JO 



of the Lyapunov exponent is the same for almost all the initial or final conditions x [|32| , |33 . 
We will show later, however, that the deviations that could occur in sets of zero measure 
may have some observable consequences. 

Taking the limit <5x — > and with substitution of fll4] ) in Eq. ([13]), one obtains 

6C{r, t; <5x) w <Jx • |* V5[r(t)]e (A - 6)(f - t) rft . (15) 

Rewriting <5x = n|5x| so that n is a unit vector, one finds the directional derivative along 
the direction of n as 

rt 
10 

If A < b this derivative remains finite in the t — > oo limit and the asymptotic field Coo(x) = 
C(x, t — > oo) is smooth ( different iable). Otherwise the derivatives of C diverge as ~ e^ x ^ 1 
leading to a nowhere-differentiable irregular asymptotic field. The exception occurs when 
<5x points along the expanding direction of the forward flow: A should be replaced by —A 
and the directional derivative is always finite. The irregular field, thus, has a filamental 
character. It should be noted that the limiting distribution is not a steady field, but 
one following the time dependence of the flow. For time-periodic flows v(x, t), will also 
be time periodic. Its singular characteristics however do not change in time. 

The characterization of the singular asymptotic field has been performed in 0] . Defining 
SCoo = <5C(x, i — > oo) and taking into account that the exponential separation (|14}) can only 



persist until it reaches the characteristic scale of the flow field, one finds (|0| in the \Sx\ — > 
limit the scaling 

6C DB (T',6x)~\6x\ a , (17) 
with a Holder exponent a given by 

a = min|pl|. (18) 

This general relation expresses the local Holder exponent in terms of the local infinite-time 
Lyapunov exponent and the chemical decay rate. As stated before, for hyperbolic systems 



this Lyapunov exponent has the value Ao everywhere but in a set of zero measure. By 
decreasing b or increasing A one finds a transition from a smooth a = 1 distribution to a 
rough a < 1 one with filamental structure. This is the so-called smooth-filamental transition 
f39|j . For uniform A = Ao, the relationship between the Holder exponent and the power 
spectrum decay exponent {(3 in the power spectrum scaling S(k) ~ k' 13 ) is /3 = 1 + 2a. 
Thus the result ([18]) implies an exponent (3 larger than one, as it is usually observed. The 
Batchelor result |35| S{k) ~ k , valid for the nondecaying scalar in smooth velocity fields 
is recovered for 6 = 0. 




FIG. 3. Snapshots of chemical patterns of the forced decaying scalar under the flow (JlQj) - Darker grey 
levels indicate smaller concentrations. The lower panels are horizontal cuts along the line y = 0.25. Left: 
b = 4.0 and T — 1.0, so that Ao = 2.67 < b. A smooth pattern is seen, in agreement with the theoretical 
arguments. Right: b = 0.1 and T = 1.0 so that Xq > b and a filamental pattern is obtained. The horizontal 
cut clearly displays the fractal nature of the field. 

In Fig. H we present snapshots of the asymptotic field evolving according to (H) and 
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(p|). For the flow we take a simple time-periodic velocity field defined in the unit square 
with periodic boundary conditions by 

v x (x, y, t) = — Y®(~2 ~ 1 m ° d T ) cos ( 27r ^) 

v y (x, y, t) = —-^-<d(t mod T — —J cos(27rx) (19) 

where Q(x) is the Heavyside step function. In our simulations U = 1.2, which produces 
a flow with a single connected chaotic region in the advection dynamics. The value of 
the numerically obtained Lyapunov exponent is Ao ~ 2.67/T. We use the source term 
S(x, y) = 1 + sin (2ttx) sin (27cy). Backward trajectories with initial coordinates on a square 
grid were calculated and used to obtain the chemical field at each point by integrating ( |iT|) 
forward in time. The smooth and the filamental behavior are found by changing the decay 
rate b, in complete agreement with the theory. 



B. Smooth- filamental transition in a nonlinear plankton dynamics model 

As we mentioned before, all the results presented in the previous section are valid for 
a more complex but stable chemical or biological dynamics, as far as the decay rate b is 
replaced the chemical Lyapunov exponent |Ac|- I n this subsection we illustrate the smooth- 
filamental transition with a simple model of plankton dynamics immersed in a meandering 
jet flow ||44|| . Additional results for nonlinear models can be found in [13. 



As in the previous subsection, diffusion is neglected so that a Lagrangian representation 
in terms of a chemical and a transport dynamical system is possible. The chemical plankton 
model is the one used in [3D|. It is similar to (|T|), but with phytoplankton growth ruled by 



a logistic term (^), and a relaxational dynamics imposed on the carrying capacity C: 

An 

§ = «*-«•. (20) 
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In the absence of flow, the only stable fixed point of model p0| ) is given by C* = Co(x) 
P* = C 5/(5 + C ), and Z* = P*/5. 
The transport dynamical system is 



dx dip 

dt dy 

dy dip 

dt dx 



(21) 



given in terms of the following streamfunction |45 
tp(x, y) — 1 — tanh 



^ y — B{t) cos [k(x — ct)\ ^ 



(22) 

I (1 + k 2 B(t) 2 sin 2 [k(x - ct)}) 2 I 
It describes a jet flowing eastwards, with meanders in the North-South direction. These 
meanders are also advected by the jet at a phase velocity c. B(t) and k are the (properly 
adimensionalized) amplitude and wavenumber of the undulation in the streamfunction. As 
we are interested in a closed flow, we impose periodic boundary conditions at the ends of 
the interval — L x < x < L x . Particles leaving the region through the right boundary are 
reinjected from the left. Nutrients are injected in and out continuously from fluid elements 
as these traverse the different regions of the carrying capacity source 

C (x, y) = l + A sm(2irx/L x ) sin(27ry/L y ). (23) 

Chaotic advection appears in this model if B is made to vary in time, for example 
periodically: B(t) = B + ecos(u;t + 6). In our calculations we use the parameter values 
B = 1.2, c = 0.12, k = 2tc/L x , L x = 7.5, L y = 4.0, u = 0.4, e = 0.3, and 9 = f . 
These values guarantee the existence of 'large scale chaos', i.e, the possibility that a test 
particle crosses the jet passing from North to South or viceversa. This is weaker than the 
requirement of hyperbolicity, but is enough to illustrate the general aspects of our theory. In 
the chemical subsystem we use A = 0.2, and 5 = 2.0, and vary the value of 7, which controls 



the relaxational time-scale. It turns out that the resulting chemical dynamics is stable 

Fig. [| shows snapshots of the long-time phytoplankton distributions. The smooth- 
filamental transition is clearly observed. More quantitative discussions can be found in 
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reference 



PI 




FIG. 4. Snapshots of phytoplankton patterns from ( pfj| ) and The lower panels are transects along 

y = 0.8. Left: 7 = 0.25; a smooth pattern is obtained. Right: for smaller relaxation rate 7 = 0.025, a fractal 
filamental pattern is generated. 



C. Open flows 

A finite region is said to be traversed by an open flow if almost all fluid elements enter 
the region and leave it forever in a finite time. A prototype is a stream passing around a 
cylindrical body. If the inflow velocity is high enough vortices form in the wake of the cylinder 
and make the flow time- dependent in this region, while the flow remains steady in front of 
the cylinder or in the far downstream region. Advected particles enter the unsteady region, 



undergo transient chaotic motion f4q- fl9f| , and finally escape and move away downstream on 
simple orbits. The time spent in the mixing region, however, depends strongly on the initial 
coordinates, with singularities on a fractal set corresponding to particles trapped forever 
in the mixing region. This is due to the existence of a non-attracting chaotic saddle (a 
fractal object of zero measure) formed by an infinite number of bounded hyperbolic orbits 
in the mixing region. The stable manifold of this chaotic saddle contains orbits coming 
from the inflow region but never escaping from the mixing zone. Thus, permanent chaotic 
advection is restricted to this fractal set of zero Lebesgue measure. Points close to the 
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unstable manifold of the chaotic saddle have spent a long time in the mixing region of the 
flow moving near chaotic orbits with a positive Lyapunov exponent. For points precisely 
on this unstable manifold, the backwards trajectories (the ones from which the Lyapunov 



exponent in (|T8| ) should be computed) remain in the chaotic saddle, thus leading to Ao > 0. 
The other trajectories spend in the mixing region just a finite time (both in the forward 
as in the backwards time direction), so that they can not contribute to the development of 
singularities in chemical distributions ((|l6l) is singular only in the i — > oo limit). In fact 
almost all trajectories are characterized by a long-time Lyapunov exponent equal to zero. 

Thus open flows provide a rather clear example of strong space-dependence of Lyapunov 
exponents. According to Eq. fll8|), the Holder exponent may be different from 1 only on the 
unstable manifold of the chaotic saddle, thus implying that the transition from smooth to 
filamental structure now only takes place in this fractal set of zero measure. The background 
chemical field is always smooth, independently of the value of b. 



Patterns of chemically decaying substances in a cylinder wake are presented in [3D]. Here 
we present the case of the stable biological model (£20J) under the jet flow (]2"T|)-(|2"2l). It is made 
open simply by not imposing the periodic boundary conditions of the previous subsection. 
The chaotic mixing region is not restricted to a finite region, but we restrict the source 
forcing ( [23]) to — L x < x < L x . This is the region shown in Fig. |], being Co = outside. We 
let the fluid particles to enter this region with very small C, P and Z concentrations. Most of 
them remain in the forcing region only for a finite time, so that, as in the cylinder wake case, 
they do not develop singularities. Only particles remaining in the chaotic saddle, i.e. the set 
of nonscaping orbits, forever will lead to diverging gradients. Fig. [| shows a phytoplankton 
pattern for the same parameter values as before and a = 0.025. Smooth zones coexist with 
singular features. The behavior is distinct from the closed flow case (Fig. ^) and represents 
a new scenario for the development of chemical or biological patterns. 
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FIG. 5. Phytoplankton filamental pattern in the open jet, with a transect along y = 0.8. 

Since the irregularities now appear only on a set of measure zero, one could ask if 
they can have any significant effect on measurable quantities. In order to clarify this, we 
consider for simplicity the case of a single chemical reactive C, and instead of the previous 
characterization of the point- wise strength of the singularities by the Holder exponent, let us 
investigate the scaling with distance 5x = \5~x\ of the spatial average of the differences SCoo. 
For simplicity, let us assume that, on the saddle, there is no distribution in the local infinite- 
time Lyapunov exponents, i.e., that the advection on the chaotic saddle is characterized by 
a single Lyapunov exponent Aq. In this case the partial fractal dimension (i.e. the dimension 



of intersections of the set with a line) of the manifolds of the chaotic saddle is [ |32| , [47|1 

/) I - — . (24) 



K 

V 



Here k is the escape rate, that is the rate of the exponential decay (~ e~ Kt ) of the number 
of fluid elements spending time longer than t in the mixing or in the forcing region. On a 
one-dimensional transect of unit length the total number of segments of length Sx is (Sx)^ 1 
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while the number of segments containing parts of the unstable manifold (with partial fractal 
dimension D) is ~ 5x~ D . Thus, according to ( [Op the spatial average of SC^ along this line, 
(<5Coo(x; 5x)), can be written as 

(<50o(x; Sx)) = {5x){5x)- £) (5x) b/Xo + 

+(Sx)[(5x)- 1 - (5x)- 3 ](5x) (25) 

where the first term pertains to the singular component (a = fe/Ao), while the second one 
pertains to the smooth component (a = 1). In the limit 5x — > the dominating behavior is 

(<fCoo(x; Sx)) ~ Sx ( , (26) 

with 

C = min{l,l + 6/A - J D} = min|l,^^| (27) 

showing that if D < (or, b + k > X ) the average will be dominated by the smooth 
component, but if the fractal dimension of the singular set is large enough it contributes 
to the scaling of (SC^r; Sr)) . We see that moments of SC may be sensible to fractal 
inhomogeneity in a, or intermittency. 

D. Anomalous scaling of structure functions 

The strongly intermittent structure of singularities in open flows is an extreme example. 
There are additional inhomogeneities affecting both to the open and to the closed flows: 
although, in the long-time limit the Lyapunov exponent is the same for almost all trajectories 
in an ergodic region, deviations can persist on fractal sets of measure zero, and as we 
saw above such sets can contribute significantly to the global scaling. The origin of these 
inhomogeneities can be traced back by analyzing the finite-time distribution of Lyapunov 
exponents. In general, the finite-time stretching rates, or local Lyapunov exponents |[33|| , 



have a certain distribution around the most probable value. This distribution approaches 
the time-asymptotic form [ 32|J33 : 



25 



P(A,t) ~t 1/2 e- G(A)t (28) 

where G(\) is a function characteristic to the advection dynamics, with the property that 
G(\q) = G'(Xq) = and G(X) > 0, and A is the most probable value of the Lyapunov 
exponent. At infinitely- long times all the measure becomes concentrated at this single value 
Ao, as stated before. There are however fractal sets of zero measure that do not share this 



unique value. The partial dimension of the set with Lyapunov exponent value A is |40 
D(\) = 1 — G(X)/X. The coexistence of such different and interwoven fractal sets is the 
signature of multifractality. 

For a robust quantitative characterization of the filamental structures, accessible to mea- 
surements and sensible to the intermittency features, we consider now the scaling properties 
of the structure functions associated with the chemical fields. For a single species C, the 
gth order structure function is defined as 

S q {Sx) = (\6C 0O (^6x)\") (29) 

where (...) represents averaging over different locations x, and q is a parameter (we will only 
consider structure functions of positive order (q > 0)). In the absence of any characteristic 
length over a certain range of scales the structure functions are expected to exhibit, as 
5x — > 0, a power-law dependence 

S q {5x) ~ 5x^ (30) 

characterized by the set of scaling exponents ( q . The scaling exponent of the power spectrum 
is given by f3 = 1 + £2- 

If the Holder exponent of the field has the same value everywhere, given by (|T^) with 
A = Ao, the scaling exponents of the resulting mono-affine field are simply 

b 

Cg = q®o = q-r-- (31) 

(we have assumed b < Ao). In general, the fractal sets of partial dimensions -D(A) should be 
taken into account in an average such as ( f?5D but for a continuous range of values of A. The 
result for the scaling exponents in (|30|) is |40[| : 
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C q = mm i q, 1 + — - D(X) S = mm ^ g, ^ 



(32) 



Equation fl27|) is a particular case of (|3^) for g = 1 and in the approximation of considering 
a single value of A on the chaotic saddle. 
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FIG. 6. The scaling exponents £ g for the structure functions of the linearly decaying scalar under the flow 
(|l9|). Same parameters as in Fig. || except b = 1.0. Thick line: the mono-fractal approximation C, q = qb/Xo. 
Thin lines: the curves Cq — Q an d Cq — [lb + C(A)]/A, for different values of A; for the numerical values 
of G(A) see |40|| . According to Eq. (|3^), the actual values of the scaling exponents are given by the lower 
envelope of this set of curves. This is confirmed by the numerically determined values of ( q (crosses). Dashed 
line: the approximation (|33[). 

Multifractality thus originates anomalous scaling (nonlinear dependence of ( q on q) of the 
structure functions. Multifractality and anomalous scaling have been observed in real data 
from plankton and temperature distributions (see for example ]51]|JFT| J5[] . Here we present 



the scaling exponents for structure functions of the linearly decaying scalar (|TT1) under the 
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flow (|19|). Fig. | shows the numerically calculated scaling exponents (i.e. obtained by direct 
application of Eq. (p9|)), and the family of lines corresponding to (^) for different values of 
A. Eq. (|32| ) predicts that the actual values of ( q are given by the lower envelope of the set 



of lines, a fact confirmed by the numerical results. The values of G(X) can be found in |10 

Also shown in the same figure are the mono-fractal approximation (( q = qb/X ), which 
appears to be accurate for small q, and the function 



c 9 

that results from a parabolic (7(A): 



■fnii \ . 

(33) 



'A \ , 2qb A 



\\A) ~A "A 



G(A) = ^£ . (34) 

This can be thought as the first term in a Taylor expansion around Ao- It gives a good 
approximation to obtain the small-g scaling exponents. Expression ( |3"3"D becomes exact for 
the Kraichnan flow 1 14111. 



V. EXCITABLE DYNAMICS 



Excitable dynamics |2C] refers to the class of dynamical systems in which one can identify 
activator and inhibitor variables with the following properties: The activator displays some 
kind of autocatalitic growth behavior, but the presence of the inhibitor controls it so that 
the dynamical system has a stable fixed point as unique global attractor. The essence of the 
excitability phenomenon is the presence of a threshold, such that if there is a perturbation 
above it the system variables reach the stable fixed point only after a large excursion in 
phase space. This behavior usually appears when the activator has a temporal response 
much faster than the inhibitor, which then takes some time before stopping the growth of 
the activator. 

Many chemical systems behave in this way, including the famous Belousov-Zhabotinsky 
reaction in adequate concentration ranges, or the electrochemical reactions occurring in 
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nerves and muscles [pO| . Truscott and Brindley |52]]T8|1 identified phytoplankton as the fast 



activator and zooplankton as the slow inhibitor in models of the type (|1|), and demonstrated 
excitable behavior when Hollings-III grazing functions are used. This opens the possibility 
of understanding some plankton dynamics features in terms of well established scenarios 



found for excitable chemical reactions. In particular explained features of red tides in 
terms of the essentials of excitable dynamics. In addition, transport processes coupled to 
excitable dynamics lead to a rich variety of pattern forming phenomena. The most widely 
studied is the appearance of excitable waves, of linear, circular, or spiral shape, in excitable 
media with diffusion [E01. 



In the next subsection we summarize the results of p7j concerning the generation of 
plankton patchiness (via plankton blooms) in excitable media without transport processes, 
and in subsection [VB| we will report some recent results for excitable dynamics under chaotic 



advection |53| , p5| , |54 



A. Localized blooms from excitable dynamics 

Excitability is a threshold phenomenon: whenever it is crossed, the system undergoes 
a large excursion in state space, after which it returns back to the unexcited state. In an 
extended system, different regions may be in different phases of the excitation cycle. If 
transport processes are not efficient enough to couple different regions, excitation will occur 
only in the places where the excitation threshold has been reached, whereas the rest of 
the system will remain essentially at the equilibrium concentrations. It was realized in |27j 



that this simple mechanism leads to localized blooms in plankton models. For adequate 
parameter values, the temporal evolution of the bloom has the characteristics of red tides 



We illustrate this phenomenon with the model presented in |52[ . It is of the form (|I|) but 
reduced to two species (P-Z) with phytoplankton logistic growth and Hollings-III grazing. 
In convenient adimensional units it reads 
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F P (P,Z) = aP(l-P) 



P 2 



P 2 + P 2 ' 



F Z (P,Z) =1 - 



-Z — mZ . 



(35) 



P 2 + P 2 

It presents excitable behavior in a range of parameters. In particular we show in Fig. [7] four 
stages of the evolution for a = 0.43, P = 0.053, 7 = 0.005, and m = 0.34. The system is 
distributed in the one-dimensional line, but each point evolves with the same dynamics (|35|) 
independently of the others. Its evolution thus depends only on its own initial conditions. 
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FIG. 7. Four stages of a plankton bloom obtained from model (35) starting from the initial condition in 
the first panel. 



The first panel in Fig. [7] shows this initial condition: Phytoplankton is everywhere at 
the equilibrium fixed-point concentration, but zooplankton presents values somehow smaller 
than equilibrium in the central zone. In the region where this depletion is larger than a 
threshold value, fast phytoplankton growth occurs, followed by growth in zooplankton and 
subsequent recovery of the equilibrium values. The bloom is localized since there is neither 
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advection nor diffusion to communicate excitation to neighboring places. Related behavior 
occurs if the initial perturbation consists in a local phytoplankton increase, or in a local 
increase of the phytoplankton growth rate a. 

B. Excitable dynamics under chaotic stirring 

Except in the most quiescent situations, fluid flow will disturb the localized blooms of 
the previous section. Chaotic advection will deform the static excited patches into long 
and thin filaments subjected to stretching and folding. Reference |53| studied this process 
including also a small diffusion. The chemical dynamics considered was the FitzHugh- 
Nagumo system, a model originally developed in the context on nerve excitation, but it is 
clear that the different scenarios found would apply to any excitable chemical or biological 
dynamics in a chaotic flow. 

Three qualitatively different kinds of behavior appeared. In the case in which the stirring 
time-scale (measured for example by the flow Lyapunov exponent) is much slower than 
the time scale for the excitation-deexcitation cycle, a localized perturbation propagates 
essentially as a deformed version of the circular waves that would propagate diffusively in a 
quiescent medium. In the other extreme case, when the stirring time-scale is faster than the 
time scale for growth of the activator, stretching of the filaments and the associated lateral 
contraction occur too fast. The largest gradients associated to filament narrowing increase 
the diffusive flux out of them, and they can not be filled-up with new excitation fast enough. 
As a result the stretched filaments become increasingly diluted and concentrations may fall 
below the excitation threshold. At this point the excitation will disappear. This means 
that a too vigorous stirring will eliminate a small localized excitation. Perhaps the most 
surprising behavior appears in a range of intermediate stirring speeds, such as the dilution 
effect stops the growth of the slow inhibitor, but not of the fast activator: The activator 
concentration inside the filament remains saturated, the lateral thinning becomes stopped 
at a finite width by the combined effect of excitation and diffusion, and the filament length 
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increases continuously as a consequence of chaotic stretching. After some time, in a closed 
flow after repeated folding, the filament will fill-up the whole system. The effect of chaotic 
advection has changed the character of an initial excitation from localized to global. After 
the whole system is excited, deexcitation proceeds as in a homogeneous well-mixed medium. 
The process is illustrated in Fig. |8], which is obtained |55|] by integration of advection- 



react ion- diffusion equations (§), with the reaction terms as in the previous subsection (Eq. 
(j35j)). Advection and diffusion for a nutrient field is also included (with Fn = 0). The 
coupling of the nutrient with phytoplankton is made via a dependence of the growth rate a 
on the nutrient concentration (a = ao(l + N/Nq)). For the flow (closed and incompressible) 
a set of randomly seeded eddies has been used |30|] . The inset in the second panel of 



Fig. |8| shows a chlorophyll filament observed from the satellite sensor SeaWiFS 42 days after 



a fertilization experiment in the Antartic [25|. Reference |55j interprets this filament as 



the result of stretching and folding of the initial excitation introduced in the fertilization 
experiment. This would support the interpretation of the observed dynamics in terms of 
the concepts of excitability and chaotic advection. The phenomenon reported here, i.e. the 
occurrence of a global excitation by the effect of chaotic advection on a localized patch, 
gives a warning about the possibility of responses of unexpected extension arising from very 
localized perturbations. 

The essential object in this kind of phenomena is the stretched filament. In |]53| , |54 



simplified equations describing the transverse filament profile are discussed. They are es- 
sentially of the form (|9|) adapted to the multicomponent case, and with the simple linear 
chemistry term n(t)P replaced by the excitable chemistry. These filament equations can be 
used to obtain analytical understanding of the transitions between the different regimes, and 



to get quantitative results. For example, the width of the excited filament is given |53|J54 



by w = cy/Da/X, where c is a numerical constant depending on the particular excitable 
model used, D is the diffusion coefficient, a the phytoplankton linear growth rate, and A the 
exponential contraction given by the strain, that can be identified with the flow Lyapunov 
exponent. 
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FIG. 8. Snapshots of phytoplankton concentration (darker grey levels indicate smaller concentrations) 
from an excitable plankton model. Time runs from top to bottom, and then from left to right. Periodic 
boundary conditions are imposed to have a closed flow. The initial localized excited patch is stretched and 
folded until the full system becomes excited. After this (not shown) homogeneous deexcitation occurs. The 
inset in the second panel shows the chlorophyll filament observed from a satellite sensor (SeaWiFS) 42 days 
after a fertilization experiment in the Antartic j25| that can be considered as the seeding of a localized 
perturbation. 
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VI. SOME WARNINGS FROM INDIVIDUAL MODELING 



In the present Lectures we have restricted our considerations to continuous differential 
equation models describing continuous concentration variables. The use of maps 0] is an 
alternative to the use of continuous differential equations. But there is an aspect of real 
chemical or biological interacting systems that is missed in both approaches: the discrete 
nature of the interacting entities (chemical molecules or the individual organisms). On gen- 
eral grounds, one would expect just a more noisy dynamics in models in which individuals 
are resolved. There are cases however in which individual quantization produces more pro- 
found effects. We summarize here the results of two recent references [j57j|58j that highlight 
this fact in very simple examples. 

Reference |57j discusses the following chemical scheme: 



A + C -► 2C 

C -> (36) 

The first autocatalitic equation can be interpreted as a simple process in which the C particle 
eats 'food' A to self-replicate (the amount of A is kept constant), and the second represents 
the death of C . The continuous description in terms of a reaction-diffusion equation for the 
concentration C would be 

= (XA - y)C + DV 2 C (37) 

A and \i are the rates of the first and second reaction in (]36f) respectively. The continuous 
model predicts extinction of C if the effective growth rate XA — \i is negative. 

Alternatively one can model the process by a large number of random walkers, each 
representing a molecule of A or C, that undergo the transformations fl3"6|) following some 
Monte Carlo dynamics. It happens p7[] that if the effective growth rate is not too negative, 



the amount of C increases instead of decreasing. The consideration of the discrete nature 
of the individuals has changed death to life. What happens is that the C individuals die in 
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most of space, but they survive if close enough to 'food' particles. The effect is a diverging 
clustering of C particles around A random walkers. This strong inhomogeneity is missed 
in the continuous model, that gives the wrong prediction. The root of the problem is in 
the use of mean multiplicative growth rates in (|37|), whereas in fact they vary randomly in 
space and time, reflecting the positions of the A particles. It is known that mild randomness 
in the multiplicative growth factors in linear equations such as (|37|) give rise to strongly 
intermittent distributions |]59| , |60| . Intermittency reflects here in the very inhomogeneous 
clustering of the surviving particles. 

A related phenomenon is observed in |58]. There, the chemical scheme is 



C ^2C 

C -> (38) 



which leads to essentially the same continuous description (p7|) , with \A replaced by A. 
Here again growth is observed in discrete simulations when the continuous model indicates 
extinction. Here there are no 'food' particles to cluster around, but clustering occurs around 
the initially seeded C particles. There are position correlations between a parent particle 
and its descendants (the ones appearing via the first autocatalitic reaction in (|3"8"|)), that are 
completely neglected in the continuous model. 



The authors of |58[ check that the reproductive correlations persist in the presence of 
chaotic advection: the clusters become filaments, but the anomalous survival remains. The 
authors are able to describe the correct behavior via continuous equations, but not for the 
average concentration, but for binary correlations. Whether or not a given discrete process 
can be modeled in a continuum model, and up to which level of statistics one should retain 
randomness and correlations, is a major question in reacting particle models and in general 
ecosystem dynamics. 
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VII. CONCLUSIONS 



We have presented a number of processes, inspired by concepts in Nonlinear Dynamics 
such as chaos and excitability, that can be useful to understand generic behaviors in chemical 
or biological systems in fluid flows. 

For the case of the linearly decaying tracer, a rather complete description is available. 
The scaling exponents of the set of structure functions, which display anomalous (multi- 
fractal) behavior, can be expressed in terms of the decay-time constant and the probability 
distribution of the Lagrangian finite-time Lyapunov exponents of the flow. The decay-time 
constant makes the tracer power spectrum steeper than the Batchelor law obtained for 
the passive tracer. The differences between open and closed flows are important, and the 
relevance of these results for nonlinear plankton models has been be pointed out. 

The study of excitable-type reactions is motivated by the presence of excitability in 
some plankton dynamics models. It contains the essentials of the observed characteristics 
of plankton blooms. Important changes in the response of excitable models are found when 
they are placed in a chaotic fluid flow. This may be relevant for the outcome of ocean 
fertilization experiments. 

Finally, some warnings have been given about the difficulties in modeling discrete in- 
dividuals (such as planktonic organisms) in terms of continuous concentration fields. Ag- 
gregation behavior, reproductive correlations, and other phenomena turn out to be difficult 
to introduce in standard continuous models, whereas they are quite naturally present in 
individual-based descriptions. 
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